{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c106181e",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "# A short tutorial on USTC-Pickers (<3 min)\n",
    "* If any question, please feel free to contact me via email: `zhujun2316@mail.ustc.edu.cn`\n",
    "* Dependencies: Obspy, PyTorch, [SeisBench v0.5.0](https://github.com/seisbench/seisbench/releases/tag/v0.5.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb1098a5",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Load a model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4606d474",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.0\n",
      "  warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plase specify a region or province to pick phases, e.g. Beijing. 42 pickers are available in the subfolder \"USTC-Pickers/model_list/v0.1\"China\n",
      "You are using the China picker, \"diting\", published by SeisBench\n"
     ]
    }
   ],
   "source": [
    "import seisbench\n",
    "import seisbench.models as sbm\n",
    "import os\n",
    "import glob\n",
    "import torch\n",
    "device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')\n",
    "from config import (en2cn, model_list)\n",
    "sample_rate = 50\n",
    "# 为拾取的目标区域选定合适的picker，如Sichuan（四川）、CSES（实验场）、TP（青藏高原）、China（中国）\n",
    "location = input('Plase specify a region or province to pick phases, e.g. Beijing. 42 pickers are available in the subfolder \"%s\"'%os.path.join('USTC-Pickers', 'model_list', 'v0.1'))\n",
    "if location not in en2cn:\n",
    "    exit('The region you specified is not available. Plase choose a region or province from below--------\\n%s\\n-----------------------------------------------------------------------------------------------'%(', '.join(en2cn)))\n",
    "elif location!='China':\n",
    "    model_save_path = glob.glob(os.path.join(model_list, '*'+en2cn[location]+'.pt'))[0]\n",
    "    print('You are using the picker located at %s\\n'%model_save_path)\n",
    "    # 模型初始化\n",
    "    picker = sbm.PhaseNet(sampling_rate=sample_rate)\n",
    "    # 加载模型\n",
    "    picker.load_state_dict(torch.load(model_save_path,\n",
    "\t\t\tmap_location=device).state_dict())\n",
    "else:\n",
    "    print('You are using the China picker, \"diting\", published by SeisBench')\n",
    "    picker = sbm.PhaseNet.from_pretrained('diting')\n",
    "    picker.sampling_rate=sample_rate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16005010",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Read waveforms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "261df044",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "../test_data/mseed/Beijing*.mseed\n",
      "3 Trace(s) in Stream:\n",
      "...BHZ | 2018-07-18T15:24:11.300000Z - 2018-07-18T15:27:11.280000Z | 50.0 Hz, 9000 samples\n",
      "...BHN | 2018-07-18T15:24:11.300000Z - 2018-07-18T15:27:11.280000Z | 50.0 Hz, 9000 samples\n",
      "...BHE | 2018-07-18T15:24:11.300000Z - 2018-07-18T15:27:11.280000Z | 50.0 Hz, 9000 samples\n"
     ]
    }
   ],
   "source": [
    "from obspy import read\n",
    "from config import (sac, mseed)\n",
    "\n",
    "# Read the 3-component waveforms as input\n",
    "#sac = sac.replace('Beijing', 'Sichuan')\n",
    "#mseed = mseed.replace('Beijing', 'Sichuan')\n",
    "inputfile = mseed\n",
    "print(inputfile)\n",
    "stream = read(inputfile)\n",
    "print(stream)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d40ae72",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Model response (optional)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a096abdc",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEGCAYAAABhMDI9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABOOElEQVR4nO2dd3gUVdfAfyeNhNB7J6EIUiMiKooUBQFBREWxgRUb9i4W9BVB0VdfK6IiWBELdsWG4qeIFOlFkCKR3gPpu/f74+5md7K7KWRDkuX8nmefnbll5uzszJx7z7n3XDHGoCiKoiheospaAEVRFKV8oYpBURRFcaCKQVEURXGgikFRFEVxoIpBURRFcRBT1gKUlDp16pikpKSyFkNRFKVCsXDhwl3GmLrB8iq8YkhKSmLBggVlLYaiKEqFQkQ2hcpTU5KiKIriQBWDoiiK4kAVgxIe5r4I0waXtRSKooSBCu9jUMoJs+4vawkUpdjk5OSQmppKZmZmWYtSasTHx9OkSRNiY2OLXCcsikFEpgCDgB3GmA6etFrA+0ASsBG4wBiz15N3H3AV4AJuNsbM8qQfD0wFEoCvgFuMBnNSFKWUSE1NpWrVqiQlJSEiZS1O2DHGsHv3blJTU0lOTi5yvXCZkqYC/fOl3Qv8YIxpDfzg2UdE2gHDgfaeOi+JSLSnzsvAKKC155P/mIqiKGEjMzOT2rVrR6RSABARateuXeweUVgUgzFmDrAnX/IQYJpnexpwjl/6dGNMljFmA7AO6CYiDYFqxpi5nl7Cm351FEVRSoVIVQpeDuf3labzub4xZiuA57ueJ70xsNmvXKonrbFnO396ACIySkQWiMiCnTt3hl1w5fC4v05tnv/z+bIWQ1GUElIWo5KCqS9TQHpgojGTjTFdjTFd69YNOnFPKQM+r5rI5KWTy1oMRalQREdHk5KSQufOnenSpQu//fYbABs3bqRDhw6OsmPHjuWpp57C5XKRkpLi+NSpU4cLL7wwLDKV5qik7SLS0Biz1WMm2uFJTwWa+pVrAmzxpDcJkq4oihKxJCQksHjxYgBmzZrFfffdx88//1xgnejo6Lw6AFu3bqVbt248+OCDYZGpNHsMnwEjPdsjgU/90oeLSCURScY6mf/wmJvSROQksUaxEX51FEVRIp4DBw5Qs2bNYtUxxjBy5EjuuuuugB7G4RKu4arvAb2AOiKSCjwMTABmiMhVwD/AMABjzAoRmQGsBHKBG40xLs+hrsc3XPVrz0dRFKXUeeTzFazcciCsx2zXqBoPD25fYJmMjAxSUlLIzMxk69at/Pjjj3l5f//9NykpKXn727Zt484773TUf+aZZ4iJieGmm24Km9xhUQzGmItCZJ0eovw4YFyQ9AVAeFSeoihKBcDflDR37lxGjBjB8uXLAWjZsqXDZDR27FhH3SVLlvDss88yf/78sI6u0pnPiqIoUGjL/khw8skns2vXLooy2jIjI4NLLrmEl156ifr164dVDlUMSsl4/ng4RuchKko4WL16NS6Xi9q1a5Oenl5g2TvvvJOePXsyaNCgsMuhikEpGbvXwdwXyloKRamweH0MYB3J06ZNIzo6usA6W7Zs4aWXXqJt27YOH0T79u155513SiyTKgYl/ORmQ0xcWUuhKBUCl8sVND0pKSnP1+DF38dQmmHkNOy2En7Sd5W1BIqilABVDEr4yc0qawkURSkBqhiU8PNcCvxZcjunoihlgyoGpXRY/QVkHYSMvWUtiaIoxUQVg1I6GAPPdoQnkspaEkVRiokqBqX0yPAs0TH9EnDllq0siqIUGVUMSunwl1+Yq9VfwO61ZSeLopRjDifsNsDll19O48aNycqygz127dpFUlJSWGRSxaAcGXTpbkUJijdW0pIlSxg/fjz33XdfketGR0czZcqUsMukiuFooaxNOVsXwwFdXkNRCqK4YbdvvfVWnnnmGXJzw/t868zno4El78PMUXDzn1CrRcmPZwz8/CR0Hp6XtDmm4Cn8fHI9RMfBg7oUq1JO+fpe2LYsvMds0BEGTCiwSEnCbjdr1oxTTz2Vt956i8GDB4dNbFUMRwNf3m6/t68Mj2LYu4EdvzzB6Zve5t24ODpmZ5NblJC/ruySn1tRIoyShN0GuP/++zn77LM566yzwiaTKoajgeyD4T2e283chHgAplerQsdde8J7fEUpCwpp2R8JihN220urVq1ISUlhxowZYZOj3PkYRKS/iKwRkXUicm9ZyxNZBHEAH9oNX95pA98Vld3rCjzytkIiQyqKEhz/sNvFYcyYMXmjlcJBuVIMIhINvAgMANoBF4lIu7KVKsL5dgzMfxVWflL0Ol/cWmD247Vr4i6RUIpy9OD1MaSkpHDhhRcWKex2ftq3b0+XLl3CJlN5MyV1A9YZY9YDiMh0YAh2fWglnGQfgtjK4MrxJBRjWcBChp7OTqzMuNo1eXC3LxzGmrhYKrvdND0MURUlkjncsNtTp0515H388cdhk6lc9RiAxsBmv/1UT5oDERklIgtEZEFxbHFHPcbTjt/3DzzeCP541ZdWnPViD27L2/y/yglsi45mcaVKjiIzqlV17J/fuCEDm/r9lbnZ8NsL4PZ7KF7qDn++XXQ5FEUpFcqbYgj2dgponhpjJhtjuhpjutatW/cIiBUhLHrTfu9Zb79/Gk/e5U1dADkZMLY6LJle5EPuiY7mvMYNWBpfKSDv2ZrV2R0V4hZ7rK41Y71yGmz4BbYugR0r4NMbi/GDFEUpDcqbYkgFh7WhCaCzosLFng3222s+ytjj6zHMe9nXWp95rbNeVhqs+jzkYQ9ER/NR1SoB6a/XqE6v5k0Klmn7cpg2yCoIL1/fY01diqKUCeVNMcwHWotIsojEAcOBz8pYpsjD+LmGV37q257/emDZVZ/D+Cbw/qXwt2/izb5QPYHC2DwfNv5acJl5k+C35337GXt11rSiHEHKlfPZGJMrIqOBWUA0MMUYs6KMxYocCltyc+eqwLT3L/Vt//4ytOzD9uhonqpd9Gn7DhGm9KVyUeIm+U+Ge7ot5GbC2P2HdU5FUYpHeesxYIz5yhhzjDGmpTFmXFnLE1Fk7odXeh7+DOS130LmAbYVFv6iAK5vUESf0D6/MQi5mfY7K83prFaU4rBlsZ23oxRKuVMMpcK+zXpDeNm62GmmKYgDWwPT5r5YotMvio8vWsFlM2BlPivi+CbwxsASnV85ipncE17tXdZSBFCSsNvJycl5cyC6d+8eNpkiWzG4cu0om2c7wMQgMYKO1lDQm+cVXubAFlgVxL2TuZ+p1auFX6ZgzLgsMG3z73YEFVhlP7Y6/OQJZZCTWfoyLf8Y/l1Y+ucJhTGQHiIESVaaXU5VCc2+TWUtQQAlCbs9ceJEFi9ezOLFi/MUSjiIbMUw/SLnvr8ieGMgPFLDvliKS262HdqZn/Q9sOKTI/NwblkMi9+F3KzSOf5/j4WdawLT573M94mVS+ecwZg6KDDttdNhfyo83cbu/zQeJjSHcfVh51+hj5W23f7fxf3PN/0GazwLD314BbzaxyqjUGFEvn3AniPzMH0im/+wUT6DKbrfnocnk52DBryMbwLjA6b9HD6rvoDPb4Hs9OLX3bvJXoP1P4VPnpIQSpkeKYyBLX/a61IAB/bupuaRangVQLlyPh8WW/60N2C/cdB9tE375AZY/E5g2bkv+sps8hsZM+lUuO7/in7OV06zjtqH9kCUx96esdc+sF4e2gv+I3fS90DlWs79jL3wfBcbjvqeTRBXxBduxj7bLQYbzvrBXRAdW3T5i8qCIKOUwkCaCFWL2lvb+Evw9GfaO/cz99nvhVOh9/1QyTN81hjbAAC4clbwY636Av7+AQY9Y8u/fynUbgl9H7X5bwyw3/f6+T1+Gm8/D+yAGL85HMb4THUTmsGwqdDoOKhSH2ITfOUO7bb3R0IQJ/7rfX3bd/0NiXV8+989aL9njHA6458+NvhvO1wy9sH7l9jthVOd5zIG0rba+9ZfNn/mPGm/3xwCD+6G6MN81WxbDlXq2c/hsugt+Gx0ocWe+OMJVu9ZHSLX2IZAdFzoyaDeYJVxvqHbbWu15Z5u98B+z72Tscfem5V9sZB8Ybcz2LplCz/OmJTXiA0Iu70llTuvuwwObgfgrrvu4rHHHgNsWIx33gny3gM7EtHbIBqzzXkvBiFyegzfjoHdf9vtYErBWwYcwy4B2zobW8O5mM2udb7W5U9P+NL/mecbvfNoLdsqNCZw0ftFU33bm/+wSuP/noUpA6xSeL6L/YB1Bj/e0NfK2ra84J7Auxc49/9Tx04Q27I4sGw5cNZujHG+FLonNeW3hCL6GorL7y/aVvPnt9p9r1IAmHKmb3tsdVgwxT7s719it9+9EFZ/aZci/fV/MGsMZB7w1ZkQJKDHx6Ps9/KP7X01Z6Iz/4PL4X+dYVwDcLt9ZSe2sPeM//EBln7g3J/YMvRvzQtnAqQVYTjvkukwvpm9vwrjieah8x6pYXuUE1vCx9cG5i94wzmD/T8FBITbsx5WfxU6f9Ip8FTr4D10sL9l1pjQZuGczECl8PPEEGXTwZ1rP/5DusHOq3HnQM6h4LL4RzDOPhhYP93Px7nvH0eW15S0+sfpfPP284y45SHM1iWAL+z24sWLWfznIq677DxbyTN829+U9M7rL9v3Udo222D2vyb7U33b3xRuqqr4PQZ/nu9StCGNnwZrPRh7A59yi20pvnC8L+unx+HYQVCvHUzp56w2oRn0DxKu94vboOuVdtvbAvz+Yfvt37Pw53+d7PekU+x3nwegQWd4d5ivTJcRwX0E3gli1/4CDT3H2bvJd8wyZHDTRgFp1zaox7IN/wQpHSYWvmEf9IL44jb78fLXN/bjZe4L0Ozkgo+x8pOim6aWfwidLrDmKC8Tmvru2dxs+PjqwHr/6wyjFwTO5Xi5O4yeD2u+caZ75bl9le2R/C/FEcaESafYc6YuhIQa8NWd9r5qP9T2YvM3cgA2zYXmJwe+gJdOhxNH2R5RYj04tCN4kMXcLNtI2bXG9qA+vApanW57vACXfASNUqxfq3U/qN4E3jrXV/+/7aD7TfYaHDsY6reHtd/BJ9fZ/LkvwPW/2XSw16RWC3jxhEBZZj8GPe+y22nbYIs12d2TfI6zXINO1iLgdsO2Jc68eu2cPcUtfwaep9Fx9juUUvPHc11P7tqZXXv2sXPXbqicr1W/NZ8M/o3HrIOw529n/u51UKd1YONw4Rsw+NkCxYksxQBOzejl7g2+l/HG/4MD/4au/+v/7Cc/L3e35qpgfBMiOvjeTaG72kXhx8cC07xhLULxSg/70G+eD6+fkZc8P74Sx2dmlbiLWKGipi59v+TH8JpTwsHH18C/iwLTD+0quGewdyNsX+EzH3rZ5fGnvHdh8Hr/LcC8lF+Z/f2jVQzBlALAG/1BoqBtkMVgXu0T+jxeJvWwLeVcv5fk8g992++c5yx/0yJr3vOSsQd+eMRuz381+Dle7m5Ne1/cFtpq4I8rp+Ch29uWQsPOcDDI6LwdK6FhijUrhZp8mZVmf3Owcxhj625bbnsXWxcDsHrdBlwuN7VrVid9i194+/Qgoyq90QHcLti9Nkj+QXuew/B1VXzFUK8d4Nfy/H6sb7vDefZm97ftz3/Nt92yD7Q+E765p2jn8pqiCuKeTb5ueGm31odNtaaK/OR76DsmNwPgtj17uXx/GukibIiNpVluLtXdxXvVf14l8TCFDWRzTDRNc32tmeGN6rOiUiWWbfiHHCAtKopaxZSv3DPv5cC0YErhlqXO+ye/UvDy3HHhkQsCzVgAvcfAbE+DyLgLDI0SwDH9fT2wXUEGMhSE18xaXB4roi9izwY7YbPm6QWXy99K92f/ZqjRLM/eH0CQdUvyyEqDSlXBnUNGZhYpfe0yucYYpj37SGDY7X3Be9d33XUXjz06Ns+s+MeXbxEX5+dvzNjjrHvHX0UKmCmmgg/Z7Nq1q1kw6/3gN5K/Wen1foEmmLs3WKWRdTD0aI6OF9gx9f4Mm2oVzuRezi7kHX9B1frFG/Vy/W92nkWoVl8oTn8YetxuRyZ5u+Mh8CqGzplZ1HK5mO0ZVdQ6O5uP/91WUFUA9kdF8Ud8JfqmZ/BOtSpMqF2r0DpFZdmGf5hZJZEmublc2bA+AH9s3MwDdWrxbZXEsJub0kSoYkxxgoyXDWP3W1t9MLNMi17BR/uM2WZ9GeFi1M/WvBPsfr5hnu1lzAphrz6cZ6E0SbkUFvv5PbqMhEXTWHXmDI5tXkRlUvdYiIkLrSwS61rT3a61BF0Uq1572LfR19JPrAuH8kWHjq4EVRv4htXWaA6VqsF2v7Wo6x0LOzx+zgadbM+mEFZt2sGxe3+AgU/mpYnIQmNM12DlI8P5XDtIi+vqfA7my4M4uLw9iUpV7IN4e74RCWO2w3lBuq3th9rvUT/50kb9ZB8EgCu+Dqxzx1/2HDflMyXUbw9t+sNV38HD+6xd+P4gXdfLPoHLv4SbF9v8Hp51nFMuht4PBJb3kuIzheyPispTCgBr4+LYEhPNjHwB8HZFR/FhVdszeL16VU5t3oTb69fl3xLMeA7FtuhoHqpbO08pAHRLasq3xemZBLteQVgfG0P3pKZ8XJJez43z4bzX7Yi0grh7A5zzsl0M/p6Nwcu0PjN4urdB0/mi4PkXB2nZgx1pcveGwHvscBi73yoFgNgg16teWzj5huB1e97jexaK+N8UyN0bfNu3BnGcj14AySF6VGDNS+e8CBe+bRUCwKJpvvz4mvZl3DDFmo5CERtvzWn1QqwdVrUBxCUGP0btVlapVPNrgOZXCmCvq/9ItX2bYKffe6lmklUeXvyVQv2O1q/R6DjnebwMeCIwLQSRoRiC0eR4535RhstVa2iHfoIdThbrGTkzxG+275XfOuuM3e95iPy69M2DzED0Pii1W8J1v1pn30O+hWxo2s128ao1ssNW795gnXnec7TsDUmnQq3kwGGtPe6Ae/9xHg/sbznnpbzdYKEszmzamP/UqcWN9W2oio7JzejdrAmP1KnN5phoXq3ha/H1bxqiV1UCphQyZvuQX7c3F7imQV3m+I9oGr0g4HrMrJLIZv/fOnY/JPVgQ6ztYs/J79TLT5X6gdfSS53W0PF865Ss5zdk9to5vu167WyjI+ViOww62JDUuzfAJTOgpt9AhKQezl5ubDyM/MJZLyrGvmBG5JvHcP4U+125lr3HrvUM842vbhXTPRvt/eD/u7rfBMPfDWxE3ZfPT3ddviHDY/x6mdflC4g46ic7XNhL/nu11312UMb9W62SPfsFX97g52zjyJ+rf7S/yfuc1Wjq2/Z+6rSGkSFibVap73MSHzsYTgqizGo2h6oN7fMnUfZZblCAGdjf6eylci3734A9Tv0OgFiF0Og4azYCqzjyU8uvYStRtn6s33Vz+40+S6gZ2hTk/45LDBJ6phhrrlR8H4OXB3YU3b4I0PTE4OnRsYEjmzpfDF/cDg06QDNbLzM3kxPesSMefhz2I3Ur5/sjbl1uh+I1OcHaE/1p0AHOLiQsReVacFcQh1IwoqLsCwDsg7XqM2g72DmPAsgsICJqsJflwKaNqepy2vjDaUYCeK961QLzb61fh1e32ZbV9vMm8/uix/g9IYHb9uzlyv1p9qUAcMGbMGMEiyvF8VBdOzzSYYa6/At40ppZfkysTOodK2iyaz1MGwx3rYfE2rBpLu60rfwSH0ebjB3kGWWSeljnb8Yecty5dHm7C6NTRnPt1d/Bz0/YVmjtlnD1D3by3XlB5n+M3Q9jq5MuQtxpdxPj7a3estiOegn13yT3sKN01noaJBd51spIOg2Ou9Q3LLTdOc56DTuFHqEXLL3fOJ8PrVK+/6R2S4jxKOMH8tnTG3Swv3vlp9DvP8HPd9pddhiv1zTlpe4x9tP5Isg64OvB3/A7pM6HThcGfwmH4rjL4M+3nGn5rQA1k/JVMhiCLAQTFQ21W/uculUbOvNrtXSOAqqRb3hvdKzzt/oTl+gMKx9fzfYs/c36dY7Jc0j7ZPLzHTTsXLD/w08JWHdB8YynkeFjWOAJkeBvz/S7+XNcOezM2Emj78f57IwP7LQtryDkunPJceeQEBO6ZTl56WSe/9P3cl82clnIspvTNlOvcj0q+XUBO07rWGi9cOE9V2Es2/BPnj+iPPHEjl2srBTHuSNnM+STIXnpyzb843zJrf2Ojr/dnrc749+tHJudk1em85udcfuNL18yYglR4nwhv7r0VZ778zkAPkvdQnJOrq1vDBhDx7d8ZoL8/50xhhlrZjDn3zn8r/f/iInK1+5yu+j4VgpnJp3JUz2DL9ye7comWqKJjvLr8eRk5PkOvP/P+4Pep13tdnbujnH7FGQIct25GAyx/i+X4uDKBUzYJlIezD5IlbgqhRc8HA7tso7l3vf7JqD64/ee2NDrRap26E/t2rWRYC3qjH1WOQWbEJZ10CrMAqwRO9J3cDD7IC1q+IXkcbt8JqCaLSAhhB/mwL9wcIdvv2Fn26Pw4u/frNHMMWnOy6Edq9i3/yDZUTVIbum8RwryMUROjwHgnEl2XHOXEXlJ/i/Feef/QGWvYoiJw23cvL/mfU5qeBLJ1X1d+kEzB/HvQeeQ1vwvAX+lALA3cy81430mg4XbF9K2VltOevckAE5vdjrP9n42QOSO0zqybOQydmXsompcVYfycBs3O9J3UDWuKuk56Xy+/nOGthrqOE842RxTPm+He+rZIb/T/JQC2JfkyPlPMbjlYNrUaoO71engFy7mgsYNWTbcJuzN3OtQCmAVhfd/NcbQ6U2n+eDsJo1YNtTTUhcJ5k508MCvD/DZ39ak8cFfH3BRW6ePwHge6lkbZ3FNx2toWrUplf1MBsYYjn/bmkCnnzWdpOpJJMYm2pfS2S8wq1IULLBzZi78wg5WqBVfi58v/LlAuf49+C/9P+oPBG+I/PrvrzSu0pjBnwwG4NpO1zL6uHxzfaJj2HZoG9UrVWdO6hz6NO1D7GEqiYnzJ/Lmyje5uO3F3HdicOf14JmD2XhgI79c+As14msU7wSJdeD0B4tUtMmJQ0jdtovSWCI4MzeTPZnWF7VJNtEg0X9ggKcHtn8LBa1FluOKRnIziImtDPudI7uyc4Vdmbup4XJTeX88sMORb4xhy6EtbM7YTEKdBJIpuPHgT4XvMbTo0MKsX77ekTbw44FsTtscULZ309481/46qFyLrWLo95FvstqbA96kU51OLNqxiCtnXRlQ98/L/nS0APO3wh8/9XEGt7QP1hvL3+C/C/8bcIxFly3imYXP8NbKtwLyvCwdsRQRYdXuVVzwxQVBy9x2/G2MaDeCXRm7eODXB7i207UcW+vYgBbY3C1zGfXdqJDnihSqxFZh7sVzWbl7Zd4L00tCTAJ/XPIH87fND/q/el+UuzJ20XtGYOTNL4Z+QfNq1kzwyNxH+PCvD4PW/37T99z2022OvO/O/44GiQ3YlbGLOgl1uPGHG5mTOsdRxv++enfVu4z/YzwADRIbsO3QNl7v9zrdGnZj+6HtnPHhGQSjQWIDvjv/O77Z+A0Hsw8ypOUQx0vb/1597JTHGNhioKPnEKxH+eRpTzIgeUDefo4rhy5vO0f+XdP4I3ZEz+LTfyYD8Gj3RxnaemjAsbYe3Jr3rP10/m/0+tDng5t9wWzqJDjn+ngVhxfvNXa7DS3ut4NIbuwfz/EthD7Ngs+h6PV+L05ocAITe9pZztsPbWd35m5a1WhFXE6GnZgKIc1t2w5t46FfH+LpXk9TNa5gc2cwsl3ZeUo+/+/wkuvO5bi3jguaB/Dswmd5ffnrIfP9/7dgvd83V7zJxAW+Wd75j1FQj6HCK4aE5AQz/bvpDGllW5PnfXYef+0NHUjty6FfcutPt7J2bxHt9x7ubDeF1B2JjDmrHfsy99Hj/R4BZbwXvqimm1AsuHQBXd8O+n8VSN/mfVm0fRG7M3cz+4LZQV90RyPzLp7Hie8G9ymd2/pcTm92Oit2reClJS8FLTO+x3hObXRq0P/c27ou7D9/rPt4HvgteOv42FrHMmPwjJDH+O2i3+j+XsEhlatXqs7+LPuS8ypDgJcWv8TLS5xzJ7yKBAq+V/1fJMM+HxYQRyj9nyup3GyKI+2+bvdx8bEX5+0H64kVdJ5gMnnz1+1I44z/ziG25q/EN7DzKWaePZNWNVs5yr+85GVeWmz/y7Nbnk1qWiqLdvhGal3T8RqaJTaib9PeJCTUzHuhehsHteJr5bX0/c//1PynmLZyGhN7TqR/ku2BZbmyiCIqoPcUrCEy5cwpnNDgBMe5vJzV4iwm9PBFUHAbN53f9Jktp581nfZ1fIMddmfspteMXnn7Hw7+kDa12hR4Hd8Z+A6d6vr+i1JTDCIyDBgLHAt0M8Ys8Mu7D7gKcAE3G2NmedKPB6YCCcBXwC3GGCMilYA3geOB3cCFxpiNhcmQkJxgWo21N8ar/V7lmm+vOezfczj855T/8OCvttu6bOQyhn8xnBW7ddG58sQpjU7h1y2FLCdaAs5rdSkfrXu78IIFMOu8WZz5UfDhq1d0uII3lr9RrOPVr1yfNrXaBPRQvDTMuZShJ1TLe4EG440z36BrA/veKG5jx/synb56OuPmFb7e1qPdH6V/cn/Gz5vAzHUfO/KWjFjC3sy9DJj6H3bvaE3lpEkB59qZvpPn/3yeMSeNKXaj6vNzPiepelLI3/hc7+fo3ay3I3/m2TN5bN5jLNy+kHqV6/HDMN8s7S0Ht4T8L5eNXBZSWS6+bHGeb+mbDd9w15y7HPmnNuxHn6STGHbMME5+92QO5jijOMdGxdKnWZ88/1X+3zPpjEmc0viUvP3SVAzHYqMkvALc6VUMItIOeA/oBjQCvgeOMca4ROQP4Bbgd6xieM4Y87WI3AB0MsZcJyLDgaHGmEJnffkrhlBc0/EaXl0WYhp9IZzd8uw8u3EwFl22iC5v2S720hFLC20dKUpFompcVV7o8wIjvxlZrHqv9H2F7o26l7j3DHBOq3P4ZN0nJT5OQUzsOZG7fr4rZP6Hgz/k/M/PD5n/Wr/XOLb2sbyx/A1eW+aLrvCf7o/x4G++eUbLRi7jrZVv8eT8J4MdpkhWh3cHvsvFX10cMv+JHk9wZtKZpLyVAtgeZIYnFMmAxLcZ3bs1zWpXLr0JbsaYVcaYYHPdhwDTjTFZxpgNwDqgm4g0BKoZY+Yaq5HeBM7xq+OddfIhcLoEHSZQPHLTm/PfGQXEoQnBkhFLWDZyGWcmhZiEhDVL+dtq8yuFg3+NodmBF8jee1Kxz5+fDrVTSnwMpWhkbh9I+sbQvpnZF8wusP6j3R8tMP/3i+YfllwA0/pP49rk0gmHDjDnQmcPIy07LUAp5B+tZ0zgyJ9rv7s2YACH47irxhdZptJWCkCBSgEoUCkAXP3t1Zzy3ikOpQDw6+JmHPzLF0onMzczpFIAeHTuo1z33XUFnstfKTx+6uMB+ff8cg8/p/oGJLw1wOfT/PrQpZz1RYih+n6U1gS3xoC/9zfVk9bYs50/3VHHGJML7AeCxusVkVEiskBEFgTL79O0D8tGLiNt1QQyNtlwEYc23BRU0PRN15C25iFH2it9X8mzO57W5LSQP7JZtdBDO69KegHjqsoZbRuSvTPQaXhlhysZe/JYlo5YyuBq75K2KvAPXjJiCf93wSLSVk1g7m/Wh/JUz6cCHt7cg8c49ltU9w2NO7i26KtBRQrjTi3ZUuE5e07DlRFkxT+gXuV6VIkNPczys3M+C+qA9fLVuV8x+p3QIQy+Ouf7AmXrUr8LT321k7RVEzi08Xqy9xb+kIO9lz4cFGRxHw9Tev7MR/1/o2Z8TfrWDBK80cOykcvoLM7re3D1YzzX+7mAst6RUF7i3daJP7nHdxR3XH1ZcDhOZ3+GHTOMDxZuwbh8x5m+erqjTNqqCWRt9zn5P/jrA4fZ84ehBTciejXuHzT9ltm35G3n9z0UhUIVg4h8LyLLg3yGFFQtSFqoWRZeW1ZBec5EYyYbY7oaY7q6M5uQvecUR/7/+vyPDxY4RyW5MxvjyrThnz8cbEeWHNowGld6S3BX5uDae8nacSaH1t9Moss35X3T7kOkrZpAIs4JLHd39QXeu/eEwJdvZrodaXFdz5YYVxXS1ozl/m4P8Pipj3Ntp2u5pcstnHfMeaRnu3h33j9AFGmrfROEDm24kbl/7+GeDz0vEXcCaasmcGbSmdSMr8mhjVbhdY1+ksyt5zrO/ek5n3J+q0vJ2XccJrecxKoJM2mrHwmZ5zrQhfcHvU/D+OAmxrQ1j5CbHny9gUc6f1Tgef/T/T888tlfHNp4PRlbnBFBW9dsTXL1ZDJzXBxcdzf9EicFOFabVm3K7DX2xf5S9x9JWzPWkX/GxPkYl51fU3vn04685ISTSM/2rRnStcFxZG0bwowBXzO4ySgOrr2Pg3/fESDzpDMmESVRfLEwh7RV4wMaC6fXu4Jhk+bS75k5rNxygI9/Cz5suVsNa9n9bmn+aKFCj8Y9WXzZUg6uvT+wItbMumftDaStmkC3pPrUqRJHn8qh/RvBGkr5ydx2dsi8YQ2eIW3VeM6o+oSjIZW5NbTSBjuAA8CdW5UtSwN7Ee6coq+uNuXTwPhtTy/0/afe/2rCGTeHPMaN7y4i49/hIfM7jv2WsSePLVCO295fXLCgQShUMRhjzjDGdAjyCd38sD0B/1VNmmAH66Z6tvOnO+qISAxQHSjSenxZ2weTtuYRDm0YTdqqcWTmuLjL80Lt264+GyfYUMHpG27mzd5zaFOrDWmrJuDObMLLl9g/z+TWIHt3b9xZjRjyotXY63ak0XPiTwAcH/soHwz6kLTV/+HQxut58K2azPnLjn1O2+lstV3V4Som/byRmCghIc7TzXbHc9+0KrSr1pvRx43O65G0f9hvVTETy5gOH5C2ahzuzKZc8to8vlnhDHK3LNWOPHFnNCdt1QQmX3QmGx67hHZRt+DKbEj6P5eTmePijc87kLn1QnoeE2RqfBGIOxQ4Aqe8cOjv28AEnxGbtmoct72/hIETN/DXn4FrG+QebAPu0LNp12yx0V4fH2rnlzRMdM54XbP9AO/9sRl3RnNy9ztj/V+TbEOndBv3PSanFmnp9jyv9bPmhaxdvUjL9IU36NG6Lrjj85R882rNycoVDq2/jbtSJnBOShJpf/nG4y9dNJh2D/nul7M7NwKiWLjexbvftcDkVufSLs4hkj0a9eHkRnZNied/XAcI0690mkc/+dnX4xz4XIgV84Bt//TA65NMW/U495xwD/e0mQnA5r0ZtLz/K0xu6Beny23rRkcJrepV4e8Q8RuXjlhKfGwM9WM7hDxW+sZR5Oztzufn+KK9GpedG5C9ryud6h0LCL+sjKNmfE2uaPQRaasmkLPvREa3sWH1T6sVGHzy7Ma30L/y2xxaOwZM4ATYQ+vuJW3VBFqbmzi08XqHmchR7u/b8bZ1374qeK/OZNtnc1CnRmRuHxiQ/94Z/8f8jXvJPZDC86e9xcmJdzrycw7YUUqZe7uStmqCtZCkOsPEZ2wZxsw//yV901VBZQhFaZmSPgOGi0glEUkGWgN/GGO2AmkicpLHfzAC+NSvjteYeT7woymiZ3zRg31Z/9hQ3JlNgGjaPuhbuOTVEda38txFdrzw0Jd+y3u59mtXnwEdG/LXYwNY9GBfvr/dF4jL7Tac8V9fS+O54cfRtnYbMLG4M2xrc8SUP7hjxhImfL2a3LS2AGTt6MuzH9iJJLmeB+HCrj4defrTP/PqnPVkZLtYv9M3quD723vy5pXdGH58W1Y+GniT1Eq0N+ngF/6PpHu/BKBprQRio+1f+P5lV5O+4RZch9ry+RLfhBnv7y8u+3aFeanIMOLOtnGnDq2/hZwDHcnZ598yc9q701Y5zR65afZ3Xdfe9viy9/gW4nm17xRemWPnxAw/wf5nbw5wrn/xyAxfa3njhLOYfcFs4vdfQNqaR7j+HTsk8kCmbdXffLq9D05seCLVtz9N9s7+/PnPvsDfk9GMm4+7mf90s8NKTW5NLus0kMtOag6uxLyH3v+3XXFKEn3a2hAwP63xTWx6dEgHFl22CNk4kbRV41mzdCjfLN+ed88AnNiiNgfX3U327h4MbzCJYJ31ab1/4Pj6Xcncek5e2op/DzL3b++6AFFc2u5SqifYCXpTf/UFusva6ZxbcHDt/Qyf/LsjrUOj6qzYcoAnTvXN93HtOYOp/aciIjSqkUAr1538OMzGccrv9zmpsVXK1WMa4cpoStaunhz860HS/nqQd895miEp1jqw62AWOS43z/1gh6cPSWnEtSf1Yd7F83hh0PUOM3LWjr5c/vpyPljos3bb6+6RL9MqYoBFqxvjzmjOn2POpZLU8JVf/ZhtdGb7wvOc2roO/eveHXCNvcTFRPHlSGcgzFnnzmLQ877lhnslpzDkGOciYZlbbA9uzEwbWPD9USeRm9aRjH99Y3Zy99tnw5Xemuzd1rIy7+IgC33lo0SKQUSGikgqcDLwpYjMAjDGrABmACuBb4AbjTHewPvXA69hHdJ/A95QpK8DtUVkHXA7EGL1GycdGlenVmIcUVHCJSc6bf7+L8X+7X2zDqfN3QjA3f2t7S0uJopaiXG0quezHXsn0gCserQ/0VH24dkwfiDndvEFk/tokb2JMlIv54bmn5C92xfffeTJVoFMOK8jXZv7ZiuP+2oVxz70DX2e9jmIWtWrwmme1n3luBjeuMLXGn3jihP46ubAFvwrlzpf+l/cdCpAXm+pX7v6xMUc3l/sSk92vDTLA66MJsy5cA6/3G3Hf7uzGnJDu0dYfMMUFl66kE8Gzg1Sy/cyzdw6lJx93QC45bQezLvoT7K2DyFt9aOkbx7BLVN9k52iPP93g8QGHPzrQdy5lTm0YXTe8Z73NDTqJNThoZ5X5fVCDvj1CDo09pnxUppY0+LT365x1LcIg5pfytDn7f/2wXUnIyLUrhK6ZzNm4LE0qmGdwN+v8ikGESE2KpYFY/oBwj970rnhncBoq68M70vWjrN4dfa+vLTKcb5r1aVZPab2f4M3h93KpQ2m5SnY22fY+DzPXpgCQK829p6dNndTXt3/G/UEb/admbdvcqsxb4Pt/HufiZRmNTAGctLakbb6UdJWPc5rZ4/h+Pq2x9OmflXWbE+jbuW6uP5+gpw9Pfj+xit58fQXWXTpIg5l2dfJdyu3k77xRrJ3DgCiwZVIStOajvAWf2zwGR4eHmxb2ZVjKyMibHx8GL9f/Dtfnv2j49kFeOsqe694FXPL7MDZ1DUqx7FgxC9c0/EaLmp7ERseP5tFD/blu9tOo22Dqnx6o30ZP3y6cxZ81i6rPOePsf7HVnWrkLbqcTK2DKNX3BR+WukzGc65y97vfdrWyzOF5aYn8cVopwI+sUVt+rStR+6B42iSeTu/X/w7Q1J876rnzxzLspHLHLPtQ1GiGAjGmJnAzBB544AAL6BnSGtAH9EYkwkMy59eGP5tnXFDO7Jw015qVI7lpUuOz2tlg33592lbjx9X7+BDT4sguU6gE/GlS7o4HqTLuyf5zEHYB++/F6Tw1PmdHcrjnv5tub5XS17+aV1ei/HeAcfm1fnw+u4s/3e/oxXg5ZMbTwlI692mHqse7U92rpvqle3IpyfP68TdH9mXx5KH+1E9wTmpxv9FBDgUmBfjSkCii7DUIELW9iFc2+FW3tjia4G4MhsRHV+E9YVLgatbPEvN+JrUjCfPPOgjmpZ141j+yJl08DPP/XhHT1rUXcbXy7bmtei9VI6LYcnD/ej8yLe4DrZjO8HX2f759kGcNtEZFXNwZ99ypYM6NeSm92zcmpvfC7LEI9Cijq2/xNNb7d7Sjqs4MbkW8zbsYeD/fCYc/0bE8kfOZMu+DFrVrYLbGFqN+Zpbz2hNTHSgwj/3ON//7W0Q5Licne6Vj1ozUr/2DRzp39zagzb1q7JxdzrJdXy/9eSWtTm5ZW3OaJvEeS/PZdsBuwzmOZ5zVY133oNf3dyDetXiqVetFX9e9idrtu/jrFW+Fuptfa3ZynuOW6YvBuI818I31qRNg6p8vXwbh7JySc+27snkOokkYweDTLr0eE4a/wOzPb2lc1Ia0bNNXZrWrJzXiLu8exJTf9vIJa/5zu//TvCSGJtIYs1EHh7cjkc+X2nT4qKtqc+Pt686kc+W/MuDn9p5So+d43uN3dzF5yeolRhHrcQ4vrnVN3ClSqUYLj/2OqaumkQc1RnT4zZObV2PulWt8o+KEsaf25n7Po7i8yU7+HyJ/V2XntSMZrXtizwhLpqcfSfiSk/CnV074HkH2xhesHEPJ7aw1/LZC1P4dPEW6lSpFPCfF0TEhd3+5tbTmD7q5KA3wC2nO2OFeG8gfwZ0cF68sWe3DygD9o/cOOEsFjxwBvcPbMt1Pe0olqVjz2T2nb14f9RJDoUC9sX9x/2n06Smb7jf97efRkrTGkHPkRAXnacUAC44oSlrHuvP2nEDApRCMDo1cR7XuOLzHPDBaFY1cJTVOSnN8hyoVWOrUm1vsPWyi0fWrgJi5xfA7f0KH11RpVIMax7rz+sju7Jxwlm0qGuV/4COPl/Bz3f1ytuunhDLggeco8by73sfTC9Pnu8cluzfOv1pjfU7eU1RXq44Jdmx770/p15hW6V70309Df/jVakUwzH1qxIVJcRER7FxwlnceobPJ3B7X9+2tzXspX0jn71/wrkd+fvxgVSO87UF/3uBnVl7U59WtG1QDRFxKAV/Oue7l0LRzu+cMVExtG9Yh3FD7Qv0jStOoEZl+7vzn6dr85qO5zHXo9Du/Th4kMn61ewL9evl1lHRq009hh7XhK5Jvui/Dw92rpswunfB850u755E1fgYuiXXYsWjdrTPikd8/pjqlWO59KTmvDaiK9Ou7MalJwUfwBCKO7rdyLKRy1g48v8Y3i2JJjWd99VF3QKfv0fPdrahVz56Jo8POoO142z4nfdHnUT3lrVZ/R8rb3SU5CkFsPeS9z1VHMpn1LRSonPTGlzQtQnrdx7ivVHB5xaICO9cfSL3fbyM2Xf2KvSYdapUYtRpznkSyXUSQz5g9arF83/3FGGN3BBUKmSxnLvObMPEWdZc4TU1eKmREE/qusuITvgnL5RB9t5uxNX8gwENb+CBXpdwynun5Dnx+rdvkPdifWvAWyRXT2bSj//ytjNWV7F5YcD9vLvpEAt3BI42Ht5mONUqVaNv876M/Hok6bnpnNXirLxwD0WhUkw0px9bPyB9xrUn06pelYBGQ50qlagcF016titvvyAu6No0IO3HO3o6TIP3DXT6aGomxjH0uMbM/PNffrijZ97LP3/jwes3KCo3n96aM9s3YMOuQ45GBMCXN/fg4ld/Z1CnRgwP8tI5t0sThqQ0DtpAyk9MdFRei/r7251DuF8b0ZWr31wQYMr1cnG3Zgw7vqnDrOmvoMCaz/wZ3q0pL8xel+cvK0gZg9cRT0CZZrUq88+edADu6HdMQJn85ZeNdTrmEyvFsOShfnmyiwhntAu8t8LFxglnMXHWal6c/TeDOzfKM2l6qRwX4/gvT2xRm3f9FEG4qPCxkhxhtxXAjqZqVCMh7+HzzqLsVLcTU/q9yTvz1vLMOmu1+3X4r7y2bCo3d7mBHHcO3d7pRofaKbzY+/Wgva4cl5s3nmnD83XjA/KKSkGzO/2Hd6Zlp3Eo51C+qJSlQ2aOiw8XpjKgQ4Ogtv3sXDej313EuKEd87r/+fF38AaaukLjX2/2nb1CNioijVyXmxVbDtA5SI851+Wm1RjfSoiLH+qb19vwcsK479mZZs1/oa63MYa1Ow5yTP2SzUmIRCJ/aU/FQat6VR0tsuf72BDhz/Z6lkox0VzULSkvr1qlatze9WZiomJIiElgav+pvNLvxaBKASA2Oor4MLUlutQreMH3qnFVj4hSAIiPjebSk5qHdPjGxUQxeUTXkEoBYLqnF+p1jheVJQ/Z0Sa1EuOOGqUAthcSTCl48/zJrxQAfr/POovfuyZ0ZAERUaVwGGiP4Shl+urpnNLoFJpWCzSLFMa0ia15ql7JewzpOems37+e7zZ9x6YDm+jeqDsXtAkealw5+tiRlslbczdxe99jgi+io5SIo2ehHqXIDG8bejZlYdRITCDEpPRCua6zLw5M5djKdKjTgQ51Qk9kUo5e6lWN544iDDhQwo+akpRiE3VK6Cn8hTGqU+QvHKQoFR1VDErxqVJ6ozIURSl7VDEoR5RoKXi4raIoZY8qBuWIkn9dWkVRyh/6lCrF5piaBU8UUhSlYqOKQSk2bWq14fvzC15Qxp9Lj720FKVRFCXcqGJQDov6ifUDFqEJhhjDPd3uKbScoijlB1UMSqnyxM7dhRdSFKVcoYpBKTUe3rWbAYfSy1oMRVGKiSoGRVEUxUFJV3CbKCKrRWSpiMwU8a1xJyL3icg6EVkjImf6pR8vIss8ec95lvjEswzo+570eSKSVBLZlLLjwuRB1IqvRa/0oiwIpChKeaOkPYbvgA7GmE7AX8B9ACLSDhgOtAf6Ay+J5M1sehkYhV0HurUnH+AqYK8xphXwDPBECWVTyogHThvPzxf+TB2XOy9tYs+JfDD4gzKUSlGUolIixWCM+dYY412c9HegiWd7CDDdGJNljNmAXd+5m4g0BKoZY+YaG9b1TeAcvzrTPNsfAqeLhlSMGPon9adtrbZlLYaiKEUgnD6GKwHvyhqNgc1+eametMae7fzpjjoeZbMfCLo0kYiMEpEFIrJg586dYfsBiqIoShHCbovI90Cw1VLGGGM+9ZQZA+QC73irBSlvCkgvqE5gojGTgclg12MIKbyiKIpSbApVDMaYAleRFpGRwCDgdONb9ScV8F8BpgmwxZPeJEi6f51UEYkBqgN7ivAblDLk0e6P8tBvD5W1GIqihJGSjkrqD9wDnG2M8R+w/hkw3DPSKBnrZP7DGLMVSBORkzz+gxHAp351Rnq2zwd+NBV9ebmjgKGth7Lg0gU8eNKDgZkaME9RKiQlXcHtBaAS8J3HT/y7MeY6Y8wKEZkBrMSamG40xrg8da4HpgIJWJ+E1y/xOvCWiKzD9hQOf4kx5YhSKboSNSrVCMy4Yw1kHjji8iiKUjJKpBg8Q0tD5Y0DxgVJXwAErOVojMkEhpVEHqXsiIsOXKydKvXsR1GUCoX29RVFURQHqhiUsNCieouyFkFRlDChikEJC82qNStrERRFCROqGBRFURQHqhgURVEUByUdrqooecw8eyYm+GR1RVEqEKoYlLDRqmbI0cuKolQg1JSkKIqiOFDFoCiKojiQih6OSETSgDVlLcdhUAfYVdZCHCYVVXaV+8hSUeWGiit7ceRuboypGywjEnwMa4wxXctaiOIiIgsqotxQcWVXuY8sFVVuqLiyh0tuNSUpiqIoDlQxKIqiKA4iQTFMLmsBDpOKKjdUXNlV7iNLRZUbKq7sYZG7wjufFUVRlPASCT0GRVEUJYyoYlAURVEcqGJQFEVRHKhiUBRFURyoYlAURVEcqGJQFEVRHKhiUBRFURyoYlAURVEcqGJQFEVRHKhiUBRFURyoYlAURVEchF0xiMgUEdkhIstD5IuIPCci60RkqYh08cvrLyJrPHn3hls2RVEUpXBKo8cwFehfQP4AoLXnMwp4GUBEooEXPfntgItEpF0pyKcoiqIUQNgVgzFmDrCngCJDgDeN5Xeghog0BLoB64wx640x2cB0T1lFURTlCFIWS3s2Bjb77ad60oKlnxjsACIyCtvbIDEx8fi2bduWjqSKoigRysKFC3eVpzWfJUiaKSA9MNGYyXgWpOjatatZsGBB+KRTFEU5ChCRTaHyykIxpAJN/fabAFuAuBDpiqIoyhGkLIarfgaM8IxOOgnYb4zZCswHWotIsojEAcM9ZZUS8tJP61iWur+sxVAUpYIQ9h6DiLwH9ALqiEgq8DAQC2CMmQR8BQwE1gHpwBWevFwRGQ3MAqKBKcaYFeGW72jkyW/W8CRr2DjhrLIWRVGUCkDYFYMx5qJC8g1wY4i8r7CKQ1EUpdyTk5NDamoqmZmZZS1KSOLj42nSpAmxsbFFrlMWPgZFUZSIIDU1lapVq5KUlIRIsPEzZYsxht27d5OamkpycnKR62lIjAjHdtAURSkNMjMzqV27drlUCgAiQu3atYvdo1HFEOGoXlCU0qW8KgUvhyOfKoYIx62aQVGUYqKKIcJxq15QlIgmOjqalJQUOnTowLBhw0hPTy/xMVUxRDjaY1CUyCYhIYHFixezfPly4uLimDRpUomPqaOSIhzVC4pyZHjk8xWs3HIgrMds16gaDw9uX+TyPXr0YOnSpSU+r/YYIhztMSjK0UFubi5ff/01HTt2LPGxtMcQ4ahiUJQjQ3Fa9uEkIyODlJQUwPYYrrrqqhIfUxVDhKNqQVEiG6+PIZyoKSnCMe6ylkBRlIqGKoYIR01JiqIUF1UMEY4qBkWJbA4ePBj2Y6piiHB0gpuiKMVFFUOEo0H0FEUpLqoYIhztMSiKUlzCrhhEpL+IrBGRdSJyb5D8u0RkseezXERcIlLLk7dRRJZ58haEW7ajEfUxKIpSXMI6j0FEooEXgb5AKjBfRD4zxqz0ljHGTAQmesoPBm4zxuzxO0xvY8yucMp1NKOKQVGU4hLuHkM3YJ0xZr0xJhuYDgwpoPxFwHthlkHxQ/WCoijFJdyKoTGw2W8/1ZMWgIhUBvoDH/klG+BbEVkoIqNCnURERonIAhFZsHPnzjCIHbloj0FRIptx48bRvn17OnXqREpKCvPmzSvxMcMdEiPYUkGh3kyDgV/zmZFOMcZsEZF6wHcistoYMyfggMZMBiYDdO3aVd98BaDOZ0WJXObOncsXX3zBokWLqFSpErt27SI7O7vExw23YkgFmvrtNwG2hCg7nHxmJGPMFs/3DhGZiTVNBSgGpejocFVFOUJ8fS9sWxbeYzboCAMmhMzeunUrderUoVKlSgDUqVMnLKcNtylpPtBaRJJFJA778v8sfyERqQ70BD71S0sUkarebaAfsDzM8h11aI9BUSKXfv36sXnzZo455hhuuOEGfv7557AcN6w9BmNMroiMBmYB0cAUY8wKEbnOk+9dWmgo8K0x5pBf9frATM/C1THAu8aYb8Ip39GI9hgU5QhRQMu+tKhSpQoLFy7kl19+Yfbs2Vx44YVMmDCByy+/vETHDXvYbWPMV8BX+dIm5dufCkzNl7Ye6BxueY52wt5j2LsRXj0drv4OarUI88EVRSku0dHR9OrVi169etGxY0emTZtWYsWgM58jnLCPSlr6AaTvgj/fDu9xFUUpNmvWrGHt2rV5+4sXL6Z58+YlPq4u1BPhhF0xRMfab1dOeI+rKEqxOXjwIDfddBP79u0jJiaGVq1aMXny5BIfVxVDhBN2F4NXMbhzw3xgRVGKy/HHH89vv/0W9uOqKSnCCXuPIcrTltAeg6JELKoYIpywO5+9isGtikFRIhVVDBFO6fkY1JSkKFD+h4QfjnyqGCKcsN+0UepjUBQv8fHx7N69u9wqB2MMu3fvJj4+vlj11Pkc4agpSVFKjyZNmpCamkp5DuYZHx9PkyZNilVHFUOEE/aGjHjiJKrzWVGIjY0lOTm5rMUIO2pKinDC7mNwu5zfiqJEHKoYIpywKwbjVQzaY1CUSEUVQ4QTdlOScdtvNSUpSsSiiiHCKT1Tko5KUpRIRRVDhBP2UUlGFYOiRDqqGCKc8PsY1JSkKJFO2BWDiPQXkTUisk5E7g2S30tE9ovIYs/noaLWVYpP2CfeuD2KwasgFEWJOMI6j0FEooEXgb7Y9Z/ni8hnxpiV+Yr+YowZdJh1lWLgDvf722tKonzO9FQUpeSEu8fQDVhnjFlvjMkGpgNDjkBdJQSl5nwupyEAFEUpOeFWDI2BzX77qZ60/JwsIktE5GsRaV/MuojIKBFZICILyvNU9PJA+J3PakpSlEgn3IpBgqTlfzUtApobYzoDzwOfFKOuTTRmsjGmqzGma926dQ9X1qOCsPsYjM54VpRIJ9yKIRVo6rffBNjiX8AYc8AYc9Cz/RUQKyJ1ilJXKT5hN/i41cegKJFOuBXDfKC1iCSLSBwwHPjMv4CINBCxkdhEpJtHht1FqasUn/APV/UcT/WCokQsYR2VZIzJFZHRwCwgGphijFkhItd58icB5wPXi0gukAEMN9beEbRuOOU7Gim1CW7qY1CUiCXsYbc95qGv8qVN8tt+AXihqHWVkhH+eQyqGBQl0tGZzxFOqc181pAYihKxqGKIcEptgpsqBkWJWFQxRDi6UI+iKMVFFUOEU2rrMWiPQVEiFlUMEY76GBRFKS6qGCKcsA9X1YV6FCXiUcUQ4ZTems/qY1CUSEUVQ4QT/lhJakpSlEhHFUOEU2qxklQxKErEooohwnGH28mQZ0rSpT0VJVJRxRDhhD9WkjeInlsX61GUCEUVQ4RTahPcQM1JihKhqGKIcMI/wU0Vg6JEOqoYIhyXn2YIywgl/x6DS/0MihKJqGKIcNwOxRCGA/qH29Yeg6JEJKoYIhz/UUlhsSqpKUlRIp6wKwYR6S8ia0RknYjcGyT/EhFZ6vn8JiKd/fI2isgyEVksIgvCLdvRiMu/gR9uU5IqBkWJSMK6gpuIRAMvAn2BVGC+iHxmjFnpV2wD0NMYs1dEBgCTgRP98nsbY3aFU66jGX8fQ9gVg/oYFCUiCXePoRuwzhiz3hiTDUwHhvgXMMb8ZozZ69n9HWgSZhkUPxympLD4GLTHoCiRTrgVQ2Ngs99+qictFFcBX/vtG+BbEVkoIqNCVRKRUSKyQEQW7Ny5s0QCRzphdz6rKUlRIp6wmpIACZIW9HUkIr2xiuFUv+RTjDFbRKQe8J2IrDbGzAk4oDGTsSYounbtqtNvCyDspiTtMShKxBPuHkMq0NRvvwmwJX8hEekEvAYMMcbs9qYbY7Z4vncAM7GmKaUE+JuSwuNj8PNmq49BUSKScCuG+UBrEUkWkThgOPCZfwERaQZ8DFxmjPnLLz1RRKp6t4F+wPIwy3fU4RyVFIYDOnoMuiaDokQiYTUlGWNyRWQ0MAuIBqYYY1aIyHWe/EnAQ0Bt4CURAcg1xnQF6gMzPWkxwLvGmG/CKd/RiKOXEHYfg/YYFCUSCbePAWPMV8BX+dIm+W1fDVwdpN56oHP+dKVkuMM+XDUXYuIhN1N9DIoSoejM5wjHFW4fg3FBdCXPwbXHoCiRiCqGCMfZYwjHAd0Q41EM6mNQlIhEFUOE43LESgpTjyFPMWiPQVEiEVUMEY7/qKSwTXDLUwzqY1CUSEQVQ4RjSmOCm/oYFCWiUcUQ4bjC7mNwqY9BUSIcVQwRjsPHEK4eg5qSFCWiUcUQ4ZRKED11PitKRKOKIcIJ+zwGt8tOcAPtMShKhKKKIcIplVhJ3h6DSxWDokQiqhgiHGPC7GNw+41K0h6DokQkqhgiHP9RSa5wdBlys6BSFbutPgZFiUhUMUQ4/sogK9ddQMkiYAzkZkCcVzFoj0FRIhFVDBGOv8O5xIrBnQvG7VMM6mNQlIhEFUOE499jyC6pYsjJsN9xlUGitMegKBGKKoYIx+U2xMfavzkrt4QzlbMP2e/YyhAVoz4GRYlQwq4YRKS/iKwRkXUicm+QfBGR5zz5S0WkS1HrKsUnPdtFjYQ4IAw9hqwD9ju+ulUOWWkllE5Rwkx2Onz7APwzr6wlqdCEdQU3EYkGXgT6AqnAfBH5zBiz0q/YAKC153Mi8DJwYhHrlhxjMMbgNtb+HhsTHdbDF3ZuxzcF7TvzjHF7vk3et/GUM55yxhi7Z7zVDO6MfTSJjyb9wEHc6XsgPS6fPMWQadsyu51QA6o3hV1rIX0PVKoG0WFfDDA4xvjkLvTb7dzO3A9VG/rmYZRT7H9r70933rfdjhKIFoiNiiIqSvwr5T9KsAP7tkUAsd8iNs/twuV2kZGVDcaNGBdRgOCGqGgkOhaRaKKiooiKikaiojyHEr9jhiY310WOy012djY52Zls3LqTv/5aSd2EKDqf2Ju4uFjioiEuShDc5Lpyyc1x4XLlkpudiduViys3B1d2OtWis0lwZ0D2QSTrAFHZB2HXX/DvQtiyCBa+CSM/hb2bYOkM2L4M2gyE0+7GxFTCbTzXWcTeUhiM8ZPf81uiowTxu3a+3xqVt20Q3MYQJUJUlOCp4Xc9Ao9b2LUqayQsY9u9BxM5GRhrjDnTs38fgDFmvF+ZV4CfjDHvefbXAL2ApMLqBqNr165mwYIFIfPnTR/PiasnhMx3GyGXKNxE4fJ0oHx/p8n7Dpbm+O350gWIkvBd2/JClomlj3mJu+QtzpE5AOSYaHI8bYzD+cX2FWXyrrO9hr7tKEzYr2WuifL7V4uPV5pgRwh1bzjTAomU+8XlecEe6WcgzSSwn0S+dJ3ExdE/UFWsT2y7qUEuMTSWXUdMluLgNr47xHtP+r69SEBa/rIAS5qN4OSrnirSeUVkoTGma7C8cDfzGgOb/fZTsb2Cwso0LmJdAERkFDDKs5slIstLIHNZUQcon3dqoVxYZ2jFlL2iXnOVu0h4TJ2s57qg6cWigl7zp+tw9dNFlbt5qIxwK4ZgDaH8TYZQZYpS1yYaMxmYDCAiC0JpvfJMRZUbKq7sKveRpaLKDRVX9nDJHW7FkAo09dtvAmwpYpm4ItRVFEVRSplwj0qaD7QWkWQRiQOGA5/lK/MZMMIzOukkYL8xZmsR6yqKoiilTFh7DMaYXBEZDcwCooEpxpgVInKdJ38S8BUwEFgHpANXFFS3CKedHM7fcASpqHJDxZVd5T6yVFS5oeLKHha5wzoqSVEURan46MxnRVEUxYEqBkVRFMVBhVYMFSWEhog0FZHZIrJKRFaIyC2e9LEi8q+ILPZ8Bpa1rPkRkY0isswj3wJPWi0R+U5E1nq+a5a1nP6ISBu/a7pYRA6IyK3l9XqLyBQR2eE/H6egaywi93nu+TUicmbZSB1S7okistoT7mamiNTwpCeJSIbftZ9UzuQOeW+U8+v9vp/MG0VksSe9ZNfbG1qhon2wDuq/gRbYoa5LgHZlLVcIWRsCXTzbVYG/gHbAWODOspavENk3AnXypT0J3OvZvhd4oqzlLOQ+2YadzFMurzdwGtAFWF7YNfbcN0uASkCy5xmILkdy9wNiPNtP+Mmd5F+uHF7voPdGeb/e+fKfBh4Kx/WuyD2GbsA6Y8x6Y0w2MB0YUsYyBcUYs9UYs8iznQasws70rqgMAaZ5tqcB55SdKIVyOvC3MWZTWQsSCmPMHGBPvuRQ13gIMN0Yk2WM2YAd3dftSMiZn2ByG2O+NcZ447H/jp2PVK4Icb1DUa6vtxexgZsuAN4Lx7kqsmIIFVqjXCMiScBxgDf842hPt3tKeTPJeDDAtyKy0BOKBKC+sXNP8HzXKzPpCmc4zoelvF9vL6GucUW6768EvvbbTxaRP0XkZxHpUVZCFUCwe6OiXO8ewHZjzFq/tMO+3hVZMRQ5hEZ5QUSqAB8BtxpjDmAjy7YEUoCt2K5geeMUY0wXbFTcG0XktLIWqKh4JkqeDXzgSaoI17swKsR9LyJjgFzgHU/SVqCZMeY44HbgXRGpVlbyBSHUvVEhrjdwEc4GUImud0VWDEUJv1FuEJFYrFJ4xxjzMYAxZrsxxmVsTO1XKaMuakEYY7Z4vncAM7EybheRhgCe7x1lJ2GBDAAWGWO2Q8W43n6Eusbl/r4XkZHAIOAS4zF4e0wxuz3bC7G2+mPKTkonBdwbFeF6xwDnAu9700p6vSuyYqgwITQ89r/XgVXGmP/6pTf0KzYUKFdRYkUkUUSqerexjsXl2Os80lNsJPBp2UhYKI5WVHm/3vkIdY0/A4aLSCURScaua/JHGcgXFBHpD9wDnG2MSfdLryt2zRVEpAVW7vVlI2UgBdwb5fp6ezgDWG2MSfUmlPh6l4V3PYxe+oHYET5/A2PKWp4C5DwV2/1cCiz2fAYCbwHLPOmfAQ3LWtZ8crfAjshYAqzwXmOgNvADsNbzXausZQ0ie2VgN1DdL61cXm+s8toK5GBbqFcVdI2BMZ57fg0woJzJvQ5rk/fe55M8Zc/z3ENLgEXA4HImd8h7ozxfb0/6VOC6fGVLdL01JIaiKIrioCKbkhRFUZRSQBWDoiiK4kAVg6IoiuJAFYOiKIriQBWDoiiK4kAVg6J4EJHaftEot/lF2zwoIi+V0jlvFZERBeQPEpFHSuPcihIKHa6qKEEQkbHAQWPMU6V4jhjsGPMuxhd4Ln8Z8ZQ5xfhNGFOU0kR7DIpSCCLSS0S+8GyPFZFpIvKtJ/79uSLypNg1K77xhD5BRI73BC9bKCKz8s2s9dIHG7Ij11PnZhFZ6QnkNh3A2JbbT9gQE4pyRFDFoCjFpyVwFjYk89vAbGNMRyADOMujHJ4HzjfGHA9MAcYFOc4pwEK//XuB44wxnYDr/NIXYKNnKsoRIaasBVCUCsjXxpgcEVmGXQjoG0/6MuwCKW2ADsB31hJENDaUQX4aYtfm8LIUeEdEPgE+8UvfATQKn/iKUjCqGBSl+GQBGGPcIpJjfI46N/aZEmCFMebkQo6TAcT77Z+FXaXrbOBBEWnvMTPFe8oqyhFBTUmKEn7WAHVF5GSwIddFpH2QcquAVp4yUUBTY8xs4G6gBlDFU+4YynckWCXCUMWgKGHG2KVmzweeEJEl2Cij3YMU/RrbQwBrbnrbY576E3jGGLPPk9cb+LI0ZVYUf3S4qqKUISIyE7jbOJdk9M+vD7xrjDn9yEqmHM2oYlCUMkRE2mDXd54TIv8EIMcYs/iICqYc1ahiUBRFURyoj0FRFEVxoIpBURRFcaCKQVEURXGgikFRFEVxoIpBURRFcfD/jDPMM/2ksEMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "response = picker.annotate(stream)\n",
    "# Plot the response\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "time_diff = response[0].stats.starttime-stream[0].stats.starttime\n",
    "fig, ax = plt.subplots(2, 1, sharex=True)\n",
    "t0 = np.arange(stream[0].stats.npts)/sample_rate\n",
    "t1 = np.arange(response[0].stats.npts)/sample_rate + time_diff\n",
    "ax[0].plot(t0, np.array([x.data for x in stream]).transpose(), label=[x.stats.channel for x in stream])\n",
    "ax[1].plot(t1, np.array([x.data for x in response[1:]]).transpose(), label=['P', 'S'])\n",
    "ax[1].set_xlabel('Time (s)')\n",
    "ax[1].set_ylim(0, 1)\n",
    "for j in ax:\n",
    "    j.legend()\n",
    "    j.set_xlim(0, t0.max()+1)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3077d818",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Retrieve the predicted picks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4e370f53",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                               fname                              p  \\\n",
      "0  ../test_data/mseed/Beijing*.mseed  [2018-07-18T15:24:51.140000Z]   \n",
      "\n",
      "                               s       p_prob       s_prob  \n",
      "0  [2018-07-18T15:25:01.800000Z]  [0.9586848]  [0.6875741]  \n"
     ]
    }
   ],
   "source": [
    "classifyoutput = picker.classify(stream, P_threshold=.3, S_threshold=.3) # Both thresholds can be tuned for your target\n",
    "picks = classifyoutput.picks\n",
    "# store the picks\n",
    "import pandas as pd\n",
    "P,S = [],[]\n",
    "for pick in picks:\n",
    "    if pick.phase=='P':\n",
    "        P.append(pick)\n",
    "    else:\n",
    "        S.append(pick)\n",
    "p, s = [str(x.peak_time) for x in P], [str(x.peak_time) for x in S]\n",
    "p_prob, s_prob = [x.peak_value for x in P], [x.peak_value for x in S]\n",
    "entry = {'fname':[inputfile], 'p':[p], 's':[s], 'p_prob':[p_prob], 's_prob':[s_prob]}\n",
    "table = pd.DataFrame(entry)\n",
    "print(table)\n",
    "table.to_csv('results/picks_demo.csv', index=False)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
